---
title: "Effects of environmental stressors on daily governance"
subtitle: "Supplementary information replication code"
author: Nick Obradovich, Dustin Tingley, and Iyad Rahwan
output: pdf_document
---

```{r opts, echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
knitr::opts_chunk$set(fig.path='.', warning=F, message=F, fig.retina=NULL, cache=T, cache.lazy = F, autodep=T, echo=F, eval=T) # set Rmarkdown options
```

```{r packages, eval=TRUE}
# packages----
library(bit64) # long integers
library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #convert to raster and vice versa
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(plyr) #for additional manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(doMC) #run multicore processes
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
library(knitr) #knitting
library(pander) #pandering
library(rasterVis) # plotting rasters
library(gridExtra) # for multiple plotting
library(GGally) # for correlation plots
library(geosphere) # spherical distance
library(Rcpp) # for C++
library(RcppArmadillo) # for C++ (and Conley errors)
library(matrixStats) # matrix calculations
library(ggExtra) # for marginal histograms
library(viridis) # for colors
library(Cairo) # pdf device
library(gtable) # for gtable filter
library(devtools)
#devtools::install_github("wilkelab/cowplot")
library(cowplot)
library(captioner)

# set cores for parallel processing----
registerDoMC(10)
```

```{r functions, eval=TRUE}
# functions for grabbing p-values and coefficients from summaries of felm models----
coef <- function(reg, name) {
    if (round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3) == 0) {
        format(round(reg$coefficients[rownames(reg$coefficients)==name, 1], 4), scientific=FALSE) 
    }
    else round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3)
}
pval <- function(reg, name) {
  if (round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3) >=0.001) 
  {round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3)} 
  else paste("< 0.001")
}

# captioner setting----
fig_num <- captioner(prefix="Fig.")
tab_num <- captioner(prefix="Table")
eq_num <- captioner(prefix="Equation")
sfig_num <- captioner(prefix="Fig. S", auto_space=FALSE)
stab_num <- captioner(prefix="Table S", auto_space=FALSE)
seq_num <- captioner(prefix="Equation S", auto_space=FALSE)
```

```{r data}
# public data----
load("./replication.RData")

# proprietary data----
# d <- readRDS('./hazel_weather_inspection_balanced_panel.rds.gz')
# hazel <- readRDS('./hazel_geolocated_master.rds.gz')
# hazel_county <- readRDS('./hazel_counties.rds.gz')
```

```{r models}
# police stops----
pol <- summary(readRDS("./log_num_stop_felm.rds.gz"))

# fatal crashes----
cr <- summary(readRDS("./prob_crash_felm.rds.gz"))

# food safety inspections----
insp <- readRDS("./facility_day_felm_quadratic_summary.rds.gz")

# food safety violations----
viol <- summary(readRDS("./num_violation_felm_quad.rds.gz"))
```

<!--
Rscript -e "rmarkdown::render('./si_replication_code.Rmd')"
-->

<!-- figure labels-->

```{r}
decon <- sfig_num(name = "decon", caption = "")
```

```{r}
margn <- sfig_num(name = "margn", caption = "")
```

```{r}
bins <- sfig_num(name = "bins", caption = "")
```

```{r}
target <- sfig_num(name = "target", caption = "")
```

```{r}
alc <- sfig_num(name = "alc", caption = "")
```

```{r}
trim.marg <- sfig_num(name = "trim.marg", caption = "")
```

```{r}
annual4.5 <- sfig_num(name = "annual4.5", caption = "")
```

```{r}
pol.grids4.5 <- sfig_num(name = "pol.grids4.5", caption = "")
```

```{r}
cr.grids4.5 <- sfig_num(name = "cr.grids4.5", caption = "")
```

```{r}
insp.grids4.5 <- sfig_num(name = "insp.grids4.5", caption = "")
```

```{r}
viol.grids4.5 <- sfig_num(name = "viol.grids4.5", caption = "")
```

```{r}
pol.grids8.5 <- sfig_num(name = "pol.grids8.5", caption = "")
```

```{r}
cr.grids8.5 <- sfig_num(name = "cr.grids8.5", caption = "")
```

```{r}
insp.grids8.5 <- sfig_num(name = "insp.grids8.5", caption = "")
```

```{r}
viol.grids8.5 <- sfig_num(name = "viol.grids8.5", caption = "")
```

<!-- table labels-->

```{r}
pol.reg <- stab_num(name = "pol.reg", caption = "")
```

```{r}
cr.reg <- stab_num(name = "cr.reg", caption = "")
```

```{r}
insp.reg <- stab_num(name = "insp.reg", caption = "")
```

```{r}
county.insp.reg <- stab_num(name = "county.insp.reg", caption = "")
```

```{r}
viol.reg <- stab_num(name = "viol.reg", caption = "")
```

```{r}
pr.viol.reg <- stab_num(name = "pr.viol.reg", caption = "")
```

```{r}
county.viol.reg <- stab_num(name = "county.viol.reg", caption = "")
```

```{r}
trimmed.tab <- stab_num(name = "trimmed.tab", caption = "")
```

# Deconvolved temperature

```{r fig.decon, echo=TRUE, eval=FALSE}
# example unit to demonstrate effect of deconvolving fixed effects on temperature data----
de <- as.data.table(demeanlist(mtx = list(crash$tmax), fl = list(crash$county_fips, crash$styrmon, as.factor(crash$date))))
setnames(de, c('tmax.demean.all'))
demeaned <- cbind(crash, de)
de <- as.data.table(demeanlist(mtx = list(crash$tmax), fl = list(crash$county_fips)))
setnames(de, c('tmax.demean.unitonly'))
demeaned <- cbind(demeaned, de)

# Washoe County NV----
p1 <- ggplot(demeaned[county_fips==32031,], aes(x=date)) + 
    ylab("Max. Temp. in \u2103") +
    xlab("Washoe County, NV") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=tmax, color='black'), size=0.15, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.15, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color='green'), size=0.15, alpha=1) +
    scale_color_manual(limits=c("black","blue4","green"), values=c("black","blue4","green"), labels=c("No Fixed Effects", "County Fixed Effects", "County, Date, and State*Month Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.position = c(0.05,0.15),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )

# St. Louis County MN----
p2 <- ggplot(demeaned[county_fips=="27139",], aes(x=date)) + 
    ylab("Max. Temp. in \u2103") +
    xlab("St. Louis County, MN") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=tmax, color='black'), size=0.15, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.15, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color='green'), size=0.15, alpha=1) +
    scale_color_manual(limits=c("black","blue4","green"), values=c("black","blue4","green"), labels=c("No Fixed Effects", "County Fixed Effects", "County, Date, and State*Month Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.position = "none",
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
# Steele County, ND----
p3 <- ggplot(demeaned[county_fips=="38093",], aes(x=date)) + 
    ylab("Max. Temp. in \u2103") +
    xlab("Steele County, ND") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=tmax, color='black'), size=0.15, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.15, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color='green'), size=0.15, alpha=1) +
    scale_color_manual(limits=c("black","blue4","green"), values=c("black","blue4","green"), labels=c("No Fixed Effects", "County Fixed Effects", "County, Date, and State*Month Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.position = "none",
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )

# Maricopa County, AZ----
p4 <- ggplot(demeaned[county_fips=="04013",], aes(x=date)) + 
    ylab("Max. Temp. in \u2103") +
    xlab("Maricopa County, AZ") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=tmax, color='black'), size=0.15, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.15, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color='green'), size=0.15, alpha=1) +
    scale_color_manual(limits=c("black","blue4","green"), values=c("black","blue4","green"), labels=c("No Fixed Effects", "County Fixed Effects", "County, Date, and State*Month Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.position = "none",
          legend.title = element_blank(),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )

# justy----
test <- ggplot(demeaned[county_fips=="04013",], aes(x=date)) + 
    ylab("Max. Temp. in \u2103") +
    xlab("Maricopa County, AZ") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=tmax, color='black'), size=0.15, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.15, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color='green'), size=0.15, alpha=1) +
    scale_color_manual(limits=c("black","blue4","green"), values=c("black","blue4","green"), labels=c("No Fixed Effects", "County Fixed Effects", "County, Date, and State*Month Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_text(size=25, angle=90, face="bold"),
          axis.text.y = element_text(size=20, face="bold"),
          axis.text.x = element_text(size=20, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position ="none",
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
justy <- gtable_filter(ggplotGrob(test), 'axis-l|ylab', trim=F)
justy <- ggdraw(justy)

# plot----
jpeg(filename = "./deconvolved_tmax.jpeg", width=18, height=10, units = "in", type="cairo", res=300)
grid.arrange(
    grobs = list(
        arrangeGrob(grobs=list(justy, arrangeGrob(p1, p2, ncol=2)), widths=c(0.05, 0.95), ncol=2, nrow=1),
        arrangeGrob(grobs=list(justy, arrangeGrob(p3, p4, ncol=2)), widths=c(0.05, 0.95), ncol=2, nrow=1)
    )
)
grid.text("a", x=unit(0.065, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.55, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("c", x=unit(0.065, "npc"), y=unit(0.47, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("d", x=unit(0.55, "npc"), y=unit(0.47, "npc"), gp=gpar(fontsize=33,font=2))
dev.off()
```

```{r decons, out.width="100%", fig.cap=paste0(sfig_num("decon", display="cite"),". Deconvolved daily temperature variation. Each panel of this figure depicts three series. The first is the raw daily temperature observations, the second is the remaining temperature variation after deconvolving the raw data around county-level fixed effects, the third is the remaining variation after deconvolving the raw data around our full set of county, date, and state-by-calendar-month fixed effects. As can be seen, the remaining variation after deconvloving the fixed effects is stationary across both seasons and years. Panel (a) plots this relationship for Washoe County, NV. Panel (b) plots this relationship for St. Louis County, MN. Panel (c) plots this relationship for Steele County, ND. Panel (d) plots this relationship for Maricopa County, AZ.")}
knitr::include_graphics("deconvolved_tmax.jpeg")
```

# Marginal effects

```{r fig.marg, echo=TRUE, eval=FALSE}
# create data----

# get cluster robust variance covariance matrices
pol.vcv <- readRDS("./log_num_stop_felm.rds.gz")$clustervcv
cr.vcv <- readRDS("./prob_crash_felm.rds.gz")$clustervcv
insp.vcv <- readRDS("./facility_day_felm_quadratic_clustervcv.rds.gz")
viol.vcv <- readRDS("./num_violation_felm_quad.rds.gz")$clustervcv

# get coefficients
pol.tcoef <- pol$coef[1]
pol.t2coef <- pol$coef[2]
cr.tcoef <- cr$coef[1]
cr.t2coef <- cr$coef[2]
insp.tcoef <- insp$coef[1]
insp.t2coef <- insp$coef[2]
viol.tcoef <- viol$coef[1]
viol.t2coef <- viol$coef[2]

# calculate marginal effects
marg <- data.table(tmax = seq(from=-15, to=45, by=.001))
marg[, pol.delta := pol.tcoef + 2*pol.t2coef*tmax]
marg[, cr.delta := cr.tcoef + 2*cr.t2coef*tmax]
marg[, insp.delta := insp.tcoef + 2*insp.t2coef*tmax]
marg[, viol.delta := viol.tcoef + 2*viol.t2coef*tmax]

# compute variance of marginal effect
marg[, pol.se := sqrt(pol.vcv["tmax","tmax"] + 4*(tmax^2)*pol.vcv["tmax2","tmax2"] + 4*tmax*pol.vcv[1,2])]
marg[, cr.se := sqrt(cr.vcv["tmax","tmax"] + 4*(tmax^2)*cr.vcv["tmax2","tmax2"] + 4*tmax*cr.vcv[1,2])]
marg[, insp.se := sqrt(insp.vcv["tmax","tmax"] + 4*(tmax^2)*insp.vcv["tmax2","tmax2"] + 4*tmax*insp.vcv[1,2])]
marg[, viol.se := sqrt(viol.vcv["tmax","tmax"] + 4*(tmax^2)*viol.vcv["tmax2","tmax2"] + 4*tmax*viol.vcv[1,2])]

# 95% confidence bounds
marg[, pol.upr := pol.delta + 1.96*pol.se]
marg[, pol.lwr := pol.delta - 1.96*pol.se]
marg[, cr.upr := cr.delta + 1.96*cr.se]
marg[, cr.lwr := cr.delta - 1.96*cr.se]
marg[, insp.upr := insp.delta + 1.96*insp.se]
marg[, insp.lwr := insp.delta - 1.96*insp.se]
marg[, viol.upr := viol.delta + 1.96*viol.se]
marg[, viol.lwr := viol.delta - 1.96*viol.se]

# create plots----
pol.plot <- ggplot(data=marg) + 
                ylab("Estimated Marginal Effect of +1\u2103 Max. Temp.\n on Log Number of Police Stops") +
                xlab("Max. Temperature in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_ribbon(aes(x=tmax, ymin=pol.lwr, ymax=pol.upr), fill="cornflowerblue", alpha=0.25) + 
                geom_line(aes(x=tmax, y=pol.delta), color="cornflowerblue", size=2) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

cr.plot <- ggplot(data=marg) + 
                ylab("Estimated Marginal Effect of +1\u2103 Max. Temp.\n on Probability of Fatal Crash") +
                xlab("Max. Temperature in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_ribbon(aes(x=tmax, ymin=cr.lwr, ymax=cr.upr), fill="red4", alpha=0.25) + 
                geom_line(aes(x=tmax, y=cr.delta), color="red4", size=2) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

insp.plot <- ggplot(data=marg) + 
                ylab("Estimated Marginal Effect of +1\u2103 Max. Temp.\n on Probability of Food Safety Inspection") +
                xlab("Max. Temperature in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_ribbon(aes(x=tmax, ymin=insp.lwr, ymax=insp.upr), fill="orange", alpha=0.25) + 
                geom_line(aes(x=tmax, y=insp.delta), color="orange", size=2) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

viol.plot <- ggplot(data=marg) + 
                ylab("Estimated Marginal Effect of +1\u2103 Max. Temp.\n on Number of Violations per Inspection") +
                xlab("Max. Temperature in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_ribbon(aes(x=tmax, ymin=viol.lwr, ymax=viol.upr), fill="green4", alpha=0.25) + 
                geom_line(aes(x=tmax, y=viol.delta), color="green4", size=2) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

# plot----
jpeg(filename = "./marg_effects.jpeg", width=22, height=9, units = "in", type="cairo", res=300)
grid.arrange(pol.plot, cr.plot, insp.plot, viol.plot, ncol=4)
grid.text("a", x=unit(0.08, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.335, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("c", x=unit(0.59, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("d", x=unit(0.84, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
dev.off()
```

```{r marg, out.width="100%", fig.cap=paste0(sfig_num("margn", display="cite"),". Marginal effects and confidence intervals of quadratic estimates. This figure plots the marginal effects of the quadratic estimates from the main text. Panel (a) plots the marginal effect of a one degree increase in maximum temperature on the log number of police stops. Panel (b) plots the marginal effect of temperature on the probability of a fatal crash. Panel (c) plots the marginal effect of temperature on the probability of a food safety inspection. Panel (d) plots the marginal effect of temperature on the number of food safety violations per inspection. Shaded regions plot the 95\\% confidence intervals of the marginal effects.")}
knitr::include_graphics("./marg_effects.jpeg")
```

# Flexible functional forms

```{r fig.bins, echo=TRUE, eval=FALSE}
# load bin plot function----
source("./binned_plot_function.R")

# police----
## bin tmax
stops[, cuttmax := cut(tmax, breaks=c(-Inf,seq(-10, 40, by=5),Inf))] # create factor variable for regresison
stops[, cuttmax := relevel(cuttmax, ref="(30,35]")] # set reference category

cut.pol <- felm(log_num_stop ~ cuttmax + prcp + trange + cloud + humid + wind 
             | county_fips + date + styrmon 
             | 0 
             | state, data=stops[full==TRUE,], exactDOF=TRUE, psdef=FALSE)

pol.p <- binned.plot(felm.est = cut.pol,
             plotvar = "cuttmax",
             breaks = 5,
             omit = c(30,35),
             linecolor = "cornflowerblue",
             pointfill = "cornflowerblue",
             errorfill = "cornflowerblue",
             xlabel = "Max. Temperature in \u2103",
             ylabel = "Log(Num. Police Stops)")
hist <- axis_canvas(pol.p, axis = "x") + 
        geom_histogram(data = stops, aes(x = stops$tmax), bins=20, fill='cornflowerblue', color='gray60', alpha=0.75, size=0.1)
pol.plot <- insert_xaxis_grob(pol.p, hist, position = "bottom", height = grid::unit(0.1, "null"))
pol.plot <- ggdraw(pol.plot)

# crash----
## bin tmax
crash[, cuttmax := cut(tmax, breaks=c(-Inf,seq(-10, 40, by=5),Inf))] # create factor variable for regresison
crash[, cuttmax := relevel(cuttmax, ref="(5,10]")] # set reference category

cut.cr <- felm(any_crash ~ cuttmax + prcp + trange + cloud + humid + wind 
           | county_fips + date + styrmon 
           | 0 
           | state, data=crash, exactDOF=TRUE, psdef=FALSE)

cr.p <- binned.plot(felm.est = cut.cr,
             plotvar = "cuttmax",
             breaks = 5,
             omit = c(5,10),
             linecolor = "red4",
             pointfill = "red4",
             errorfill = "red4",
             xlabel = "Max. Temperature in \u2103",
             ylabel = "P(Fatal Crash) in Percentage Points")
hist <- axis_canvas(cr.p, axis = "x") + 
        geom_histogram(data = crash, aes(x = crash$tmax), bins=20, fill='red4', color='gray60', alpha=0.75, size=0.1)
cr.plot <- insert_xaxis_grob(cr.p, hist, position = "bottom", height = grid::unit(0.1, "null"))
cr.plot <- ggdraw(cr.plot)

# county violations----
cut.viol <- readRDS("./violation_county_binned.rds.gz")

viol.p <- binned.plot(felm.est = cut.viol,
             plotvar = "cuttmax",
             breaks = 5,
             omit = c(15,20),
             linecolor = "green4",
             pointfill = "green4",
             errorfill = "green4",
             xlabel = "Max. Temperature in \u2103",
             ylabel = "County Number of Violations per Inspection")
hist <- axis_canvas(viol.p, axis = "x") + 
        geom_histogram(data = hazel_county, aes(x = hazel_county$tmax), bins=20, fill='green4', color='gray60', alpha=0.75, size=0.1)
viol.plot <- insert_xaxis_grob(viol.p, hist, position = "bottom", height = grid::unit(0.1, "null"))
viol.plot <- ggdraw(viol.plot)

# county inspections----
cut.insp <- readRDS("./inspection_county_binned.rds.gz")

insp.p <- binned.plot(felm.est = cut.insp,
             plotvar = "cuttmax",
             breaks = 5,
             omit = c(15,20),
             linecolor = "orange",
             pointfill = "orange",
             errorfill = "orange",
             xlabel = "Max. Temperature in \u2103",
             ylabel = "County Log(Number of Inspections)")
hist <- axis_canvas(insp.p, axis = "x") + 
        geom_histogram(data = hazel_county, aes(x = hazel_county$tmax), bins=20, fill='orange', color='gray60', alpha=0.75, size=0.1)
insp.plot <- insert_xaxis_grob(insp.p, hist, position = "bottom", height = grid::unit(0.1, "null"))
insp.plot <- ggdraw(insp.plot)

# plot----
jpeg(filename = "./binned_plots.jpeg", width=20, height=9, units = "in", type="cairo", res=300)
grid.arrange(pol.plot, cr.plot, insp.plot, viol.plot, ncol=4)
grid.text("a", x=unit(0.065, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.31, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("c", x=unit(0.57, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("d", x=unit(0.815, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
dev.off()
```

```{r binned, out.width="100%", fig.cap=paste0(sfig_num("bins", display="cite"),". Binned maximum temperatures for estimating $f()$. This figure plots the marginal effects estimated from binning maximum temperatures into 5$^\\circ{}$C bins. Y-axes report the estimated change in each respective outcome variable as compared to the omitted temperature bin in each relationship. For each estimation, we omit the bin containing the extremum for that relationship. Panel (a) plots the marginal effect of maximum temperature on the log number of police stops. Panel (b) plots the marginal effect of maximum temperature on the probability of a fatal crash. Panel (c) plots the marginal effect of temperature on the county-aggregated, log number of food safety inspections. Panel (d) plots the marginal effect of temperature on the county-day average number of food safety violations per inspection. Shaded regions plot the 95\\% confidence intervals of the estimated marginal effects.")}
knitr::include_graphics("./binned_plots.jpeg")
```

# Examination of potential targeting

```{r fig.target, echo=TRUE, eval=FALSE}
# models----
q1 <- readRDS("./prob_inspection_quartile_1_viol_felm_quad_summary.rds.gz")
q2 <- readRDS("./prob_inspection_quartile_2_viol_felm_quad_summary.rds.gz")
q3 <- readRDS("./prob_inspection_quartile_3_viol_felm_quad_summary.rds.gz")
q4 <- readRDS("./prob_inspection_quartile_4_viol_felm_quad_summary.rds.gz")

# targeting plot----
t <- seq(from=-15, to=45, by=0.001)
q1hat <- q1$coef[1]*t + q1$coef[2]*t^2
q1hat <- q1hat - max(q1hat)
q2hat <- q2$coef[1]*t + q2$coef[2]*t^2
q2hat <- q2hat - max(q2hat)
q3hat <- q3$coef[1]*t + q3$coef[2]*t^2
q3hat <- q3hat - max(q3hat)
q4hat <- q4$coef[1]*t + q4$coef[2]*t^2
q4hat <- q4hat - max(q4hat)

qdata <- as.data.table(cbind(t, q1hat, q2hat, q3hat, q4hat))

q.plot <- ggplot(qdata) + 
    ylab("P(Inspection) in Percentage Points") +
    xlab("Max. Temp. in \u2103") +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=t, y=q1hat, color="springgreen3"), size=2, alpha=0.75) + 
    geom_line(aes(x=t, y=q2hat, color="darkorchid"), size=2, alpha=0.75) + 
    geom_line(aes(x=t, y=q3hat, color="goldenrod"), size=2, alpha=0.75) + 
    geom_line(aes(x=t, y=q4hat, color="gray50"), size=2, alpha=0.75) + 
    coord_cartesian(ylim=c(-0.2, 0.05), expand=TRUE) +
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    scale_color_manual(limits=c("springgreen3", "darkorchid", "goldenrod", "gray50"), values=c("springgreen3", "darkorchid", "goldenrod", "gray50"), labels=c("1st Quartile of Total Violations","2nd Quartile of Total Violations","3rd Quartile of Total Violations","4th Quartile of Total Violations")) +
    guides(color = guide_legend(override.aes = list(fill = NULL), order=1), fill = guide_legend(order=2)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, face="bold"),
          axis.title.y = element_text(size=25, face="bold"),
          axis.text.y = element_text(size=20, face="bold"),
          axis.text.x = element_text(size=20, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          legend.title=element_blank(),
          legend.margin = margin(-0.6,0,0,0, unit="cm"),
          legend.position = c(.25, 0.25),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(q.plot, axis = "x") + 
        geom_histogram(data = stops, aes(x = tmax), bins=20, fill='springgreen3', color='gray60', alpha=0.75, size=0.1)
q.plot <- insert_xaxis_grob(q.plot, hist, position = "bottom", height = grid::unit(0.1, "null"))
q.plot <- ggdraw(q.plot)

# plot----
jpeg(filename = "./targeting.jpeg", width=10, height=7, units = "in", type="cairo", res=300)
q.plot
dev.off()
```

```{r targeting, out.width="100%", fig.cap=paste0(sfig_num("target", display="cite"),". Examination of possible substitution of inspection effort from low risk to high risk facilities. This figure plots the effects of maximum temperatures on the probability of food safety inspection across quartiles of total violations by facility. This estimation follows the same procedure as in the main text Fig. 3a, but represents separately estimating the regression for each quartile of total food safety violations by facility. The first quartile represents facilities with the lowest numbers of violations in our sample while the fourth quartile represents the facilities with the highest number of violations. As can be seen, the effects of hot temperatures on probability of inspection are negative across each quartile, suggesting that there is no substantial redirection of regulatory effort from low-risk violators to high-risk violators as a function of maximum temperatures. Each relationship is significant at the p < 0.001 level.")}
knitr::include_graphics("./targeting.jpeg")
```

# Alcohol involved crashes

```{r fig.alc, echo=TRUE, eval=FALSE}
# regressions----
cr.drunk <- felm(any_drunk ~ tmax + tmax2 + prcp + trange + cloud + humid + wind 
               | county_fips + date + styrmon 
               | 0 
               | state, data=crash[(any_drunk==100 & any_crash==100) | (any_drunk==0 & any_crash==0),], exactDOF=TRUE)

cr.nodrunk <- felm(any_crash ~ tmax + tmax2 + prcp + trange + cloud + humid + wind 
                 | county_fips + date + styrmon 
                 | 0 
                 | state, data=crash[(any_drunk==0 & any_crash==100) | (any_drunk==0 & any_crash==0),], exactDOF=TRUE)

# drunk----
## convert coefficients in inspections to percentage points
t <- stops[tmax >= -15 & tmax <=45, seq(from=min(tmax, na.rm=T), to=max(tmax, na.rm=T), by=0.001)]
insp.t <- pol$coef[1]*t + pol$coef[2]*t^2
cr.t <- cr.drunk$coef[1]*t + cr.drunk$coef[2]*t^2 + 0.15 - max(insp.t) # slight re-shift for tangency
insp.t <- insp.t - max(insp.t) # put max of parabola at zero
tmax <- stops[, sample(tmax, length(t))]
t.data <- as.data.table(cbind(t, insp.t, cr.t, tmax))

drunk.plot <- ggplot(t.data) + 
    ylab("Log(Num. Police Stops) \n and P(Fatal Crash)") +
    xlab("Max. Temp. in \u2103") +
    geom_ribbon(aes(x=t, ymin=insp.t, ymax=cr.t, fill="gray85")) +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=t, y=insp.t, color="cornflowerblue"), size=3, alpha=0.75) + 
    geom_line(aes(x=t, y=cr.t, color="red4"), size=3, alpha=0.75) + 
    coord_cartesian(xlim=c(min(t), max(t)), ylim=c(-0.75, 0.8), expand=TRUE) +
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    scale_color_manual(limits=c("cornflowerblue", "red4"), values=c("cornflowerblue", "red4"), labels=c("Police Stops","P(Fatal Crash)")) +
    scale_fill_manual(limits=c("gray85"), values=c("gray85"), labels=c("Marginal Regulatory Gap")) +
    guides(color = guide_legend(override.aes = list(fill = NULL), order=1), fill = guide_legend(order=2)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.margin = margin(-0.6,0,0,0, unit="cm"),
          legend.position = c(.25, 0.75),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(drunk.plot, axis = "x") + 
        geom_histogram(data = t.data, aes(x = tmax), bins=20, fill='red4', color='gray60', alpha=0.75, size=0.1)
drunk.plot.full <- insert_xaxis_grob(drunk.plot, hist, position = "bottom", height = grid::unit(0.1, "null"))
drunk.plot.full <- ggdraw(drunk.plot.full)

# justy----
test <- ggplot(t.data) + 
    ylab("Log(Num. Police Stops) \n and P(Fatal Crash) in Percentage Points") +
    xlab("Max. Temp. in \u2103") +
    geom_ribbon(aes(x=t, ymin=insp.t, ymax=cr.t, fill="gray85")) +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=t, y=insp.t, color="cornflowerblue"), size=3, alpha=0.75) + 
    geom_line(aes(x=t, y=cr.t, color="red4"), size=3, alpha=0.75) + 
    coord_cartesian(xlim=c(min(t), max(t)), ylim=c(-0.75, 0.8), expand=TRUE) +
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    scale_color_manual(limits=c("cornflowerblue", "red4"), values=c("cornflowerblue", "red4"), labels=c("Police Stops","P(Fatal Crash)")) +
    scale_fill_manual(limits=c("gray85"), values=c("gray85"), labels=c("Marginal Regulatory Gap")) +
    guides(color = guide_legend(override.aes = list(fill = NULL), order=1), fill = guide_legend(order=2)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_text(size=22, angle=90, vjust=-0.5, face="bold"),
          axis.text.y = element_text(size=22, angle=0, face="bold"),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_line(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position = c(.75, 0.075),
          legend.text=element_text(size=16, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(test, axis = "x") + 
        geom_histogram(data = t.data, aes(x = tmax), bins=20, fill='red4', color='gray60', alpha=0.75, size=0.1)
test <- insert_xaxis_grob(test, hist, position = "bottom", height = grid::unit(0.1, "null"))
justy <- gtable_filter(test, 'axis-l|ylab', trim=F)
justy <- ggdraw(justy)

# no drunk----
## convert coefficients in inspections to percentage points
t <- stops[tmax >= -15 & tmax <=45, seq(from=min(tmax, na.rm=T), to=max(tmax, na.rm=T), by=0.001)]
insp.t <- pol$coef[1]*t + pol$coef[2]*t^2
cr.t <- cr.nodrunk$coef[1]*t + cr.nodrunk$coef[2]*t^2 + 0.17 - max(insp.t) # slight re-shift for tangency
insp.t <- insp.t - max(insp.t) # put max of parabola at zero
tmax <- stops[, sample(tmax, length(t))]
t.data <- as.data.table(cbind(t, insp.t, cr.t, tmax))

nodrunk.plot <- ggplot(t.data) + 
    ylab("Log(Num. Police Stops) \n and P(Fatal Crash)") +
    xlab("Max. Temp. in \u2103") +
    geom_ribbon(aes(x=t, ymin=insp.t, ymax=cr.t, fill="gray85")) +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=t, y=insp.t, color="cornflowerblue"), size=3, alpha=0.75) + 
    geom_line(aes(x=t, y=cr.t, color="red4"), size=3, alpha=0.75) + 
    coord_cartesian(xlim=c(min(t), max(t)), ylim=c(-0.75, 0.8), expand=TRUE) +
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    scale_color_manual(limits=c("cornflowerblue", "red4"), values=c("cornflowerblue", "red4"), labels=c("Police Stops","P(Fatal Crash)")) +
    scale_fill_manual(limits=c("gray85"), values=c("gray85"), labels=c("Marginal Regulatory Gap")) +
    guides(color = guide_legend(override.aes = list(fill = NULL), order=1), fill = guide_legend(order=2)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.margin = margin(-0.6,0,0,0, unit="cm"),
          legend.position = "none",
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(nodrunk.plot, axis = "x") + 
        geom_histogram(data = t.data, aes(x = tmax), bins=20, fill='red4', color='gray60', alpha=0.75, size=0.1)
nodrunk.plot.full <- insert_xaxis_grob(nodrunk.plot, hist, position = "bottom", height = grid::unit(0.1, "null"))
nodrunk.plot.full <- ggdraw(nodrunk.plot.full)

# plot----
jpeg(filename = "./alc_noalc.jpeg", width=15, height=7.5, units = "in", type="cairo", res=300)
grid.arrange(grobs=list(
  justy,
  arrangeGrob(drunk.plot.full, nodrunk.plot.full, ncol=2)),
  widths=c(0.08, 0.92), ncol=2) # plot all together
grid.text("a", x=unit(0.10, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.56, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("Alcohol\nInvolved\nCrashes", x=unit(0.31, "npc"), y=unit(0.32, "npc"), gp=gpar(fontsize=25,font=2))
grid.text("No Alcohol\nInvolved\nCrashes", x=unit(0.78, "npc"), y=unit(0.32, "npc"), gp=gpar(fontsize=25,font=2))
dev.off()
```

```{r alcho, out.width="100%", fig.cap=paste0(sfig_num("alc", display="cite"),". Alcohol involved versus no alcohol involved crashes as a function of maximum temperatures. This figure plots the same police stops functional form as in the main text, with separate estimates of the fatal crash and maximum temperature functional form. Panel (a) estimates the effect of maximum temperatures on the probability of fatal crashes where alcohol was involved. Panel (b) estimates the effect of maximum temperatures on the probability of fatal crashes where no alcohol was detected. The relationship for no alcohol involved crashes is steeper than for no alcohol involved crashes, though each relationship is significant at the p < 0.001 level.")}
knitr::include_graphics("./alc_noalc.jpeg")
```

# Regression tables

## Police stops

```{r pol.regs, echo=TRUE}
# balanced panel----
quad <- felm(log_num_stop ~ tmax + I(tmax^2) + prcp + trange + cloud + humid + wind 
             | county_fips + date + styrmon 
             | 0 
             | state, data=stops[full==TRUE,], exactDOF=TRUE)
heat.quad <- felm(log_num_stop ~ nwsheatindex + I(nwsheatindex^2) + prcp + trange + cloud + wind 
                  | county_fips + date + styrmon 
                  | 0 
                  | state, data=stops[full==TRUE,], exactDOF=TRUE)

# all stops----
quad.all <- felm(log_num_stop ~ tmax + I(tmax^2) + prcp + trange + cloud + humid + wind 
             | county_fips + date + styrmon 
             | 0 
             | state, data=stops, exactDOF=TRUE)
heat.quad.all <- felm(log_num_stop ~ nwsheatindex + I(nwsheatindex^2) + prcp + trange + cloud + wind 
                  | county_fips + date + styrmon 
                  | 0 
                  | state, data=stops, exactDOF=TRUE)

```

```{r pol.tab, echo=TRUE, results='asis'}
invisible(stargazer(quad, heat.quad, quad.all, heat.quad.all,
          type='latex',
          title = paste0(stab_num("pol.reg", display="cite"), '.', " Police Stops Regression Table"),
          model.numbers = TRUE,
          column.labels = c('Restricted', 'Restricted','Full', 'Full'),
          dep.var.labels.include = FALSE,
          dep.var.caption = 'Dependent Variable: Log Number of Stops',
          covariate.labels = c("TMAX", "TMAX²", "HEAT INDEX", "HEAT INDEX²","PRCP", "TRANGE", "CLOUD", "HUMID", "WIND"),
          add.lines = list(c("County FE", "Yes", "Yes", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes", "Yes", "Yes"), 
                           c("State:Month FE", "Yes", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          star.cutoffs=c(0.01, 0.005, 0.001)
          ))
```

## Fatal crashes

```{r cr.regs, echo=TRUE}
cr <- felm(any_crash ~ tmax + tmax2 + prcp + trange + cloud + humid + wind 
           | county_fips + date + styrmon 
           | 0 
           | state, data=crash, exactDOF=TRUE)
cr.heat <- felm(any_crash ~ nwsheatindex + nwsheatindex2 + prcp + trange + cloud + wind 
                | county_fips + date + styrmon 
                | 0 
                | state, data=crash, exactDOF=TRUE)
cr.drunk <- felm(any_drunk ~ tmax + tmax2 + prcp + trange + cloud + humid + wind 
               | county_fips + date + styrmon 
               | 0 
               | state, data=crash[(any_drunk==100 & any_crash==100) | (any_drunk==0 & any_crash==0),], exactDOF=TRUE)

cr.nodrunk <- felm(any_crash ~ tmax + tmax2 + prcp + trange + cloud + humid + wind 
                 | county_fips + date + styrmon 
                 | 0 
                 | state, data=crash[(any_drunk==0 & any_crash==100) | (any_drunk==0 & any_crash==0),], exactDOF=TRUE)
```

```{r cr.tab, echo=TRUE, results='asis'}
invisible(stargazer(cr, cr.heat, cr.drunk, cr.nodrunk,
          type='latex',
          title = paste0(stab_num("cr.reg", display="cite"), '.', " Fatal Crashes Regression Table"),
          model.numbers = TRUE,
          column.labels = c('Full', 'Full', 'Alcohol', 'No Alcohol'),
          dep.var.labels.include = FALSE,
          dep.var.caption = 'Dependent Variable: Fatal Crash (0/1)',
          covariate.labels = c("TMAX", "TMAX²", "HEAT INDEX", "HEAT INDEX²","PRCP", "TRANGE", "CLOUD", "HUMID", "WIND"),
          add.lines = list(c("County FE", "Yes", "Yes","Yes","Yes"),
                           c("Date FE", "Yes", "Yes","Yes","Yes"), 
                           c("State:Month FE", "Yes", "Yes","Yes","Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          star.cutoffs=c(0.01, 0.005, 0.001)
          ))
```

## Food safety inspections
### Probability of food safety inspection

```{r insp, echo=TRUE, eval=FALSE, results='asis'}
insp.tmax2 <- readRDS("./facility_day_felm_quadratic.rds.gz")
insp.heat2 <- readRDS("./facility_day_felm_quadratic_heat_index.rds.gz")

invisible(stargazer(insp.tmax2, insp.heat2,
          type='latex',
          title = paste0(stab_num("insp.reg", display="cite"), '.', " Probability of Food Safety Inspection Regression Table"),
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          dep.var.caption = 'Dependent Variable: Inspection (0/1)',
          covariate.labels = c("TMAX", "TMAX²", "HEAT INDEX", "HEAT INDEX²", "PRCP", "TRANGE", "CLOUD", "HUMID", "WIND"),
          add.lines = list(c("Facility FE", "Yes","Yes"),
                           c("Date FE", "Yes","Yes"), 
                           c("State:Month FE", "Yes","Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          star.cutoffs=c(0.01, 0.005, 0.001)
          ))

rm(insp.tmax2, insp.heat2)
invisible(gc())
```

### County-aggregated inspections

```{r county.insp.regs, echo=TRUE}
num.tmax2 <- readRDS("./inspection_county_quad.rds.gz")
num.heat2 <- readRDS("./inspection_county_heat_index_quad.rds.gz")
```

```{r county.insp.tab, echo=TRUE, results='asis'}
invisible(stargazer(num.tmax2, num.heat2,
          type='latex',
          title = paste0(stab_num("county.insp.reg", display="cite"), '.', " Log Number of Food Safety Inspections in County"),
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          dep.var.caption = 'Dependent Variable: Log Number of Inspections',
          covariate.labels = c("TMAX", "TMAX²","HEAT INDEX", "HEAT INDEX²", "PRCP", "TRANGE", "CLOUD", "HUMID", "WIND"),
          add.lines = list(c("County FE", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes"), 
                           c("State:Month FE", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Regressions are weighted by number of facilities per county.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          star.cutoffs=c(0.01, 0.005, 0.001)
          ))
```

## Food safety violations
### Number of violations per inspection

```{r viol.regs, echo=TRUE}
# number of violations
num.tmax2 <- readRDS("./num_violation_felm_quad.rds.gz")
num.heat2 <- readRDS("./num_violation_felm_heat_index_quad.rds.gz")

# probability of any violation
pr.tmax2 <- readRDS("./prob_any_violation_felm_quad.rds.gz")
pr.heat2 <- readRDS("./prob_any_violation_felm_heat_index_quad.rds.gz")
```

```{r viol.tab, echo=TRUE, results='asis'}
invisible(stargazer(num.tmax2, num.heat2,
          type='latex',
          title = paste0(stab_num("viol.reg", display="cite"), '.', " Food Safety Number of Violations Regression Table"),
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          dep.var.caption = 'Dependent Variable: Number of Violations per Inspection',
          covariate.labels = c("TMAX", "TMAX²","HEAT INDEX", "HEAT INDEX²", "PRCP", "TRANGE", "CLOUD", "HUMID", "WIND"),
          add.lines = list(c("Facility FE", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes"), 
                           c("State:Month FE", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          star.cutoffs=c(0.01, 0.005, 0.001)
          ))
```

### Probability of violation per inspection

```{r pr.viol.tab, echo=TRUE, results='asis'}
invisible(stargazer(pr.tmax2, pr.heat2,
          type='latex',
          title = paste0(stab_num("pr.viol.reg", display="cite"), '.', " Food Safety Probability of Violation Regression Table"),
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          dep.var.caption = 'Dependent Variable: Any Violation (0/1)',
          covariate.labels = c("TMAX", "TMAX²","HEAT INDEX", "HEAT INDEX²", "PRCP", "TRANGE", "CLOUD", "HUMID", "WIND"),
          add.lines = list(c("Facility FE", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes"), 
                           c("State:Month FE", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          star.cutoffs=c(0.01, 0.005, 0.001)
          ))
```

### County-aggregated violations

```{r county.viol.regs, echo=TRUE}
# violations per inspection county level
num.tmax2 <- readRDS("./violation_county_quad.rds.gz")
num.heat2 <- readRDS("./violation_county_heat_index_quad.rds.gz")
```

```{r county.viol.tab, echo=TRUE, results='asis'}
invisible(stargazer(num.tmax2, num.heat2,
          type='latex',
          title = paste0(stab_num("county.viol.reg", display="cite"), '.', " Food Safety Violations per Inspection, County Mean"),
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          dep.var.caption = 'Dependent Variable: Mean Violations per Inspection',
          covariate.labels = c("TMAX", "TMAX²","HEAT INDEX", "HEAT INDEX²", "PRCP", "TRANGE", "CLOUD", "HUMID", "WIND"),
          add.lines = list(c("County FE", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes"), 
                           c("State:Month FE", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Regressions are weighted by number of inspections per county.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          star.cutoffs=c(0.01, 0.005, 0.001)
          ))
```

# Trimmed temperature

## Trimmed temperature marginal effects

```{r trim.regs, echo=TRUE}
# regressions----
pol.trim <- felm(log_num_stop ~ tmax + tmax2 + prcp + trange + cloud + humid + wind 
             | county_fips + date + styrmon 
             | 0 
             | state, data=stops[full==TRUE & tmax >= quantile(tmax, p=0.025) & tmax <= quantile(tmax, p=0.975),], exactDOF=TRUE)

cr.trim <- felm(any_crash ~ tmax + tmax2 + prcp + trange + cloud + humid + wind 
           | county_fips + date + styrmon 
           | 0 
           | state, data=crash[tmax >= quantile(tmax, p=0.025) & tmax <= quantile(tmax, p=0.975),], exactDOF=TRUE)

insp.trim <- readRDS("./facility_day_felm_quadratic_trimmed_tmax_summary.rds.gz")
viol.trim <- readRDS("./num_violation_felm_quad_trimmed_tmax.rds.gz")
```

```{r fig.trim.marg, echo=TRUE, eval=FALSE}
# create data----

# get cluster robust variance covariance matrices
pol.vcv <- pol.trim$clustervcv
cr.vcv <- cr.trim$clustervcv
insp.vcv <- readRDS("./facility_day_felm_quadratic_trimmed_tmax_clustervcv.rds.gz")
viol.vcv <- viol.trim$clustervcv

# get coefficients
pol.tcoef <- pol.trim$coef[1]
pol.t2coef <- pol.trim$coef[2]
cr.tcoef <- cr.trim$coef[1]
cr.t2coef <- cr.trim$coef[2]
insp.tcoef <- insp.trim$coef[1]
insp.t2coef <- insp.trim$coef[2]
viol.tcoef <- viol.trim$coef[1]
viol.t2coef <- viol.trim$coef[2]

# calculate marginal effects
marg <- data.table(tmax = seq(from=-15, to=45, by=.001))
marg[, pol.delta := pol.tcoef + 2*pol.t2coef*tmax]
marg[, cr.delta := cr.tcoef + 2*cr.t2coef*tmax]
marg[, insp.delta := insp.tcoef + 2*insp.t2coef*tmax]
marg[, viol.delta := viol.tcoef + 2*viol.t2coef*tmax]

# compute variance of marginal effect
marg[, pol.se := sqrt(pol.vcv["tmax","tmax"] + 4*(tmax^2)*pol.vcv["tmax2","tmax2"] + 4*tmax*pol.vcv[1,2])]
marg[, cr.se := sqrt(cr.vcv["tmax","tmax"] + 4*(tmax^2)*cr.vcv["tmax2","tmax2"] + 4*tmax*cr.vcv[1,2])]
marg[, insp.se := sqrt(insp.vcv["tmax","tmax"] + 4*(tmax^2)*insp.vcv["tmax2","tmax2"] + 4*tmax*insp.vcv[1,2])]
marg[, viol.se := sqrt(viol.vcv["tmax","tmax"] + 4*(tmax^2)*viol.vcv["tmax2","tmax2"] + 4*tmax*viol.vcv[1,2])]

# 95% confidence bounds
marg[, pol.upr := pol.delta + 1.96*pol.se]
marg[, pol.lwr := pol.delta - 1.96*pol.se]
marg[, cr.upr := cr.delta + 1.96*cr.se]
marg[, cr.lwr := cr.delta - 1.96*cr.se]
marg[, insp.upr := insp.delta + 1.96*insp.se]
marg[, insp.lwr := insp.delta - 1.96*insp.se]
marg[, viol.upr := viol.delta + 1.96*viol.se]
marg[, viol.lwr := viol.delta - 1.96*viol.se]

# create plots----
pol.plot <- ggplot(data=marg) + 
                ylab("Estimated Marginal Effect of +1\u2103 Max. Temp.\n on Log Number of Police Stops") +
                xlab("Estimated on Trimmed\nMax. Temperature in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_ribbon(aes(x=tmax, ymin=pol.lwr, ymax=pol.upr), fill="cornflowerblue", alpha=0.25) + 
                geom_line(aes(x=tmax, y=pol.delta), color="cornflowerblue", size=2) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

cr.plot <- ggplot(data=marg) + 
                ylab("Estimated Marginal Effect of +1\u2103 Max. Temp.\n on Probability of Fatal Crash") +
                xlab("Estimated on Trimmed\nMax. Temperature in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_ribbon(aes(x=tmax, ymin=cr.lwr, ymax=cr.upr), fill="red4", alpha=0.25) + 
                geom_line(aes(x=tmax, y=cr.delta), color="red4", size=2) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

insp.plot <- ggplot(data=marg) + 
                ylab("Estimated Marginal Effect of +1\u2103 Max. Temp.\n on Probability of Food Safety Inspection") +
                xlab("Estimated on Trimmed\nMax. Temperature in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_ribbon(aes(x=tmax, ymin=insp.lwr, ymax=insp.upr), fill="orange", alpha=0.25) + 
                geom_line(aes(x=tmax, y=insp.delta), color="orange", size=2) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

viol.plot <- ggplot(data=marg) + 
                ylab("Estimated Marginal Effect of +1\u2103 Max. Temp.\n on Number of Violations per Inspection") +
                xlab("Estimated on Trimmed\nMax. Temperature in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_ribbon(aes(x=tmax, ymin=viol.lwr, ymax=viol.upr), fill="green4", alpha=0.25) + 
                geom_line(aes(x=tmax, y=viol.delta), color="green4", size=2) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

# plot----
jpeg(filename = "./marg_effects_trimmed.jpeg", width=22, height=9, units = "in", type="cairo", res=300)
grid.arrange(pol.plot, cr.plot, insp.plot, viol.plot, ncol=4)
grid.text("a", x=unit(0.08, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.335, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("c", x=unit(0.59, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("d", x=unit(0.84, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=33,font=2))
dev.off()
```

```{r trim.marg, out.width="100%", fig.cap=paste0(sfig_num("trim.marg", display="cite"),". Marginal effects and confidence intervals of quadratic estimates of trimmed maximum temperature regressions. This figure plots the marginal effects of the quadratic estimates from the trimmed maximum temperature distribution regressions. Panel (a) plots the marginal effect of a one degree increase in maximum temperature on the log number of police stops. Panel (b) plots the marginal effect of temperature on the probability of a fatal crash. Panel (c) plots the marginal effect of temperature on the probability of a food safety inspection. Panel (d) plots the marginal effect of temperature on the number of food safety violations per inspection. The marginal effects are plotted across the full distribution of temperature for comparability with the estimates from the main text regressions, even though these parameters are estimated on the constrained 2.5-97.5 percentiles of maximum temperatures. Shaded regions plot the 95\\% confidence intervals of the marginal effects.")}
knitr::include_graphics("./marg_effects_trimmed.jpeg")
```

## Trimmed temperature regression table

```{r trim.tab, echo=TRUE, eval=FALSE, results='asis'}
insp.trim <- readRDS("./facility_day_felm_quadratic_trimmed_tmax.rds.gz")

invisible(stargazer(pol.trim, cr.trim, insp.trim, viol.trim,
          type='latex',
          title = paste0(stab_num("trimmed.tab", display="cite"), '.', " Trimmed Maximum Temperature Regression Table"),
          model.numbers = TRUE,
          column.labels = c('Police Stops', 'Fatal Crashes','Food Safety\nInspections', 'Food Safety\nViolations'),
          dep.var.labels.include = FALSE,
          dep.var.caption = '',
          covariate.labels = c("TMAX", "TMAX²", "PRCP", "TRANGE", "CLOUD", "HUMID", "WIND"),
          add.lines = list(c("County FE", "Yes", "Yes", "No", "No"),
                           c("Facility FE", "No", "No", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes", "Yes", "Yes"), 
                           c("State:Month FE", "Yes", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          star.cutoffs=c(0.01, 0.005, 0.001)
          ))

rm(insp.trim)
invisible(gc())
```

# Climate impact projections

## Annualized projections, RCP4.5 emissions scenario

```{r annual.4.5, echo=TRUE, eval=FALSE}
# calculate RCP8.5 differences for plotting----
id <- c('2050','2099') # set year id

## stops
meandiff.pol <- stack(mean(100*(exp(1)^diff2050.rcp45.pol - 1)), mean(100*(exp(1)^diff2099.rcp45.pol - 1))) # change in mean daily police stops in percentage points (convert changes in fitted logs to percentage points)
meandiff.pol <- setZ(meandiff.pol, id) # set year as z variable
names(meandiff.pol) <- c('2050','2099') # set names
meandiff.pol <- trim(meandiff.pol, padding=0, values=NA) # trim outer NAs

## crashes
meandiff.cr <- stack(mean(diff2050.rcp45.cr), mean(diff2099.rcp45.cr)) # change in mean daily fatal crash risk in percentage points
meandiff.cr <- setZ(meandiff.cr, id) # set year as z variable
names(meandiff.cr) <- c('2050','2099') # set names
meandiff.cr <- trim(meandiff.cr, padding=0, values=NA) # trim outer NAs

## food inspections
meandiff.insp <- stack(mean(diff2050.rcp45.insp*1000), mean(diff2099.rcp45.insp*1000)) # change in mean daily inspections per 100k facilities (multiply percentage points by 1000)
meandiff.insp <- setZ(meandiff.insp, id) # set year as z variable
names(meandiff.insp) <- c('2050','2099') # set names
meandiff.insp <- trim(meandiff.insp, padding=0, values=NA) # trim outer NAs

##food violations
meandiff.viol <- stack(mean(diff2050.rcp45.viol*1000), mean(diff2099.rcp45.viol*1000)) # change in mean daily violations per 1k inspections
meandiff.viol <- setZ(meandiff.viol, id) # set year as z variable
names(meandiff.viol) <- c('2050','2099') # set names
meandiff.viol <- trim(meandiff.viol, padding=0, values=NA) # trim outer NAs

# police projection----
theme <- rasterTheme(region=colorRampPalette(c("orange3","orange","gray80","cornflowerblue","steelblue","blue","blue4"))(100)) # define color palette

pol.proj <- levelplot(meandiff.pol, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               names=c('2050','2099'), # set names
               layout=c(2,1), # horizontal layout
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
pol.proj <- update(pol.proj, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Police Stops Due to Future Warming\n(Percentage Points)") # remove plot borders and add subtitle below legend
pol.proj <- pol.proj + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# crash projection----
theme <- rasterTheme(region=colorRampPalette(c("gray90","lavenderblush","darksalmon","red1","red2","red3","red4"))(100)) # define color palette

cr.proj <- levelplot(meandiff.cr, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               names=c('2050','2099'), # set names
               layout=c(2,1), # horizontal layout
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
cr.proj <- update(cr.proj, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Fatal Crash Risk Due to Future Warming\n(Percentage Points)") # remove plot borders and add subtitle below legend
cr.proj <- cr.proj + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# food safety inspection projection----
theme <- rasterTheme(region=colorRampPalette(c("darkorchid4","darkorchid","gray80","yellow","yellow2","gold","gold3"))(100)) # define color palette

insp.proj <- levelplot(meandiff.insp, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               names=c('2050','2099'), # set names
               layout=c(2,1), # horizontal layout
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
insp.proj <- update(insp.proj, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Food Safety Inspections Due to Future Warming\n(per 100,000 Facilities)") # remove plot borders and add subtitle below legend
insp.proj <- insp.proj + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# food safety violation projection----
theme <- rasterTheme(region=colorRampPalette(c("gray95","darkolivegreen2","darkolivegreen3","darkolivegreen4","darkolivegreen"))(100)) # define color palette

viol.proj <- levelplot(meandiff.viol, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               names=c('2050','2099'), # set names
               layout=c(2,1), # horizontal layout
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
viol.proj <- update(viol.proj, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Food Safety Violations Due to Future Warming\n(per 1,000 Inspections)") # remove plot borders and add subtitle below legend
viol.proj <- viol.proj + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot----
jpeg(filename = './annualRCP45.jpeg', width=14, height=6.5, units = "in", type="cairo", res=300)
grid.arrange(
    arrangeGrob(pol.proj, cr.proj, ncol=2),
    rectGrob(gp=gpar(col="white")),
    arrangeGrob(insp.proj, viol.proj, ncol=2), nrow=3, heights=c(0.48, 0.04, 0.48))
grid.text("a", x=unit(0.03, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=25,font=2))
grid.text("b", x=unit(0.53, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=25,font=2))
grid.text("c", x=unit(0.03, "npc"), y=unit(0.45, "npc"), gp=gpar(fontsize=25,font=2))
grid.text("d", x=unit(0.53, "npc"), y=unit(0.45, "npc"), gp=gpar(fontsize=25,font=2))
dev.off()
```

```{r annual4.5, out.width="100%", fig.cap=paste0(sfig_num("annual4.5", display="cite"),". Projected change in regulatory behaviors and outcomes in the US due to future warming. This figure presents the 25km x 25km grid cell projections of the effects of warming in the United States over this century on the regulatory outcomes examined in this study. Projections are calculated using downscaled climatic model data from NASA's NEX under the RCP4.5 emissions scenario (RCP8.5 projections presented in main text) across the mean of the 21 CMIP5 models in the ensemble. We couple these climate model data with the estimates from our historical statistical models to project the mean effects of climate change on each outcome. Panel (a) depicts projected changes to police stops, panel (b) depicts projected changes to fatal crash risk, panel (c) depicts projected changes to food safety inspections, and panel (d) depicts projected changes to food safety violations.")}
knitr::include_graphics("./annualRCP45.jpeg")
```

## By-month grid-cell projections

### RCP4.5 emissions scenario

```{r fig.grids4.5, echo=TRUE, eval=FALSE}
# pol----
id2050 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set month id

id2099 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set month id

forecast2050 <- setZ(100*(exp(1)^diff2050.rcp45.mn.pol - 1), id2050) # set year as z variable and convert to percentage points
names(forecast2050) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set names

forecast2099 <- setZ(100*(exp(1)^diff2099.rcp45.mn.pol - 1), id2099) # set year as z variable and convert to percentage points
names(forecast2099) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set names

## trim forecast rasters to remove outer NA values
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## combine forecasts together
forecast <- stack(forecast2050, forecast2099)

## plot raster forecasts
nam <- c(paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" "), paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" "))

theme <- rasterTheme(region=colorRampPalette(c("orange3","orange","lightblue1","cornflowerblue","steelblue","blue1","blue3","blue4"))(100)) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(4,6), # horizontal layout
               names.attr=nam, # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Police Stops Due to Future Warming\n(Percentage Points)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot
jpeg(filename = "./police_forecastRCP45.jpeg", width=8, height=8, units = "in", type="cairo", res=300)
p
dev.off()

# cr----
id2050 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set month id

id2099 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set month id

forecast2050 <- setZ(diff2050.rcp45.mn.cr, id2050) # set year as z variable
names(forecast2050) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set names

forecast2099 <- setZ(diff2099.rcp45.mn.cr, id2099) # set year as z variable
names(forecast2099) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set names

## trim forecast rasters to remove outer NA valuse
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## combine forecasts together
forecast <- stack(forecast2050, forecast2099)

## plot raster forecasts
nam <- c(paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" "), paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" "))

theme <- rasterTheme(region=colorRampPalette(c("gray65","gray80","lavenderblush","darksalmon","red1","red2","red3","red4"))(100)) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(4,6), # horizontal layout
               names.attr=nam, # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Fatal Crash Risk Due to Future Warming\n(Percentage Points)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot
jpeg(filename = "./crash_forecastRCP45.jpeg", width=8, height=8, units = "in", type="cairo", res=300)
p
dev.off()

# insp----
id2050 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set month id

id2099 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set month id

forecast2050 <- setZ(diff2050.rcp45.mn.insp*1000, id2050) # set year as z variable and convert to per 100k
names(forecast2050) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set names

forecast2099 <- setZ(diff2099.rcp45.mn.insp*1000, id2099) # set year as z variable and convert to per 100k
names(forecast2099) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set names

## trim forecast rasters to remove outer NA valuse
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## combine forecasts together
forecast <- stack(forecast2050, forecast2099)

## plot raster forecasts
nam <- c(paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" "), paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" "))

theme <- rasterTheme(region=colorRampPalette(c("darkorchid4","darkorchid","gray80","yellow","yellow2","gold","gold3"))(100)) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(4,6), # horizontal layout
               names.attr=nam, # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Food Safety Inspections Due to Future Warming\n(per 100,000 Facilities)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot
jpeg(filename = "./inspection_forecastRCP45.jpeg", width=8, height=8, units = "in", type="cairo", res=300)
p
dev.off()

# viol----
id2050 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set month id

id2099 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set month id

forecast2050 <- setZ(diff2050.rcp45.mn.viol*1000, id2050) # set year as z variable and convert to per 1k
names(forecast2050) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set names

forecast2099 <- setZ(diff2099.rcp45.mn.viol*1000, id2099) # set year as z variable and convert to per 1k
names(forecast2099) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set names

## trim forecast rasters to remove outer NA valuse
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## combine forecasts together
forecast <- stack(forecast2050, forecast2099)

## plot raster forecasts
nam <- c(paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" "), paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" "))

theme <- rasterTheme(region=colorRampPalette(c("gray85","gray95","darkolivegreen1","darkolivegreen2","darkolivegreen3","darkolivegreen4","darkolivegreen"))(100)) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(4,6), # horizontal layout
               names.attr=nam, # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Food Safety Violations Due to Future Warming\n(per 1,000 Inspections)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot
jpeg(filename = "./violation_forecastRCP45.jpeg", width=8, height=8, units = "in", type="cairo", res=300)
p
dev.off()
```

```{r pol.grids4.5, out.width="100%", fig.cap=paste0(sfig_num("pol.grids4.5", display="cite"),". Grid cell projections of change in number of police stops for RCP4.5. This figure presents the 25km x 25km grid cell projections of the potential impact of warming maximum temperatures on the percentage point change in police stops in future months. In this figure, downscaled climatic model maximum temperature projections are averaged across the 21 models in the ensemble and then coupled with our historical model parameters to produce an estimated percentage point change in police stops in each geographic location for the months of 2050 and 2099. Areas of the northern United States -- where non-summer temperatures are currently coldest -- may experience increases in net police regulatory activity due to future warming. However, southern portions of the U.S. may experience net reductions in police regulatory stops due to future warming.")}
knitr::include_graphics("./police_forecastRCP45.jpeg")
```

```{r cr.grids4.5, out.width="100%", fig.cap=paste0(sfig_num("cr.grids4.5", display="cite"),". Grid cell projections of change in number of fatal vehicular crashes for RCP4.5. This figure presents the 25km x 25km grid cell projections of the potential impact of warming maximum temperatures on the percentage point change in fatal crashes in future months. In this figure, downscaled climatic model maximum temperature projections are averaged across the 21 models in the ensemble and then coupled with our historical model parameters to produce an estimated percentage point change in fatal crashes in each geographic location for the months of 2050 and 2099. Areas of the northern United States -- where non-summer temperatures are currently coldest -- may experience decreases in net fatal crashes in the winter due to future warming. However, southern portions of the U.S. -- especially in summer months -- may experience net increases in fatal vehicular crashes due to future warming.")}
knitr::include_graphics("./crash_forecastRCP45.jpeg")
```

```{r insp.grids4.5, out.width="100%", fig.cap=paste0(sfig_num("insp.grids4.5", display="cite"),". Grid cell projections of change in food safety inspections for RCP4.5. This figure presents the 25km x 25km grid cell projections of the potential impact of warming maximum temperatures on the change in food safety regulatory inspections in future months. In this figure, downscaled climatic model maximum temperature projections are averaged across the 21 models in the ensemble and then coupled with our historical model parameters to produce an estimated change in food safety inspections per 1,000 facilities in each geographic location for the months of 2050 and 2099. Areas of the northern United States may experience increases in net food safety inspection activity during the winter months due to future warming. Southern portions of the U.S. may experience net reductions in inspection activity -- particularly during summer months when food safety risks are highest -- due to future warming.")}
knitr::include_graphics("./inspection_forecastRCP45.jpeg")
```

```{r viol.grids4.5, out.width="100%", fig.cap=paste0(sfig_num("viol.grids4.5", display="cite"),". Grid cell projections of change in number of food safety violations for RCP4.5. This figure presents the 25km x 25km grid cell projections of the potential impact of warming maximum temperatures on the change in food safety violations in future months. In this figure, downscaled climatic model maximum temperature projections are averaged across the 21 models in the ensemble and then coupled with our historical model parameters to produce an estimated change in food safety violations per 1,000 inspections for each geographic location for the months of 2050 and 2099. Our projection suggests that warming may increase food safety violation rates across the entire US by 2099.")}
knitr::include_graphics("./violation_forecastRCP45.jpeg")
```

### RCP8.5 emissions scenario

```{r fig.grids8.5, echo=TRUE, eval=FALSE}
# pol----
id2050 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set month id

id2099 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set month id

forecast2050 <- setZ(100*(exp(1)^diff2050.mn.pol - 1), id2050) # set year as z variable and convert to percentage points
names(forecast2050) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set names

forecast2099 <- setZ(100*(exp(1)^diff2099.mn.pol - 1), id2099) # set year as z variable and convert to percentage points
names(forecast2099) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set names

## trim forecast rasters to remove outer NA values
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## combine forecasts together
forecast <- stack(forecast2050, forecast2099)

## plot raster forecasts
nam <- c(paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" "), paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" "))

theme <- rasterTheme(region=colorRampPalette(c("orange4","orange3","orange","lightblue1","cornflowerblue","steelblue","blue1","blue2","blue3","blue4"))(100)) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(4,6), # horizontal layout
               names.attr=nam, # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Police Stops Due to Future Warming\n(Percentage Points)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot
jpeg(filename = "./police_forecastRCP85.jpeg", width=8, height=8, units = "in", type="cairo", res=300)
p
dev.off()

# cr----
id2050 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set month id

id2099 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set month id

forecast2050 <- setZ(diff2050.mn.cr, id2050) # set year as z variable
names(forecast2050) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set names

forecast2099 <- setZ(diff2099.mn.cr, id2099) # set year as z variable
names(forecast2099) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set names

## trim forecast rasters to remove outer NA valuse
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## combine forecasts together
forecast <- stack(forecast2050, forecast2099)

## plot raster forecasts
nam <- c(paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" "), paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" "))

theme <- rasterTheme(region=colorRampPalette(c("gray65","gray80","lavenderblush","darksalmon","red1","red2","red3","red4"))(100)) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(4,6), # horizontal layout
               names.attr=nam, # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Fatal Crash Risk Due to Future Warming\n(Percentage Points)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot
jpeg(filename = "./crash_forecastRCP85.jpeg", width=8, height=8, units = "in", type="cairo", res=300)
p
dev.off()

# insp----
id2050 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set month id

id2099 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set month id

forecast2050 <- setZ(diff2050.mn.insp*1000, id2050) # set year as z variable and convert to per 100k
names(forecast2050) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set names

forecast2099 <- setZ(diff2099.mn.insp*1000, id2099) # set year as z variable and convert to per 100k
names(forecast2099) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set names

## trim forecast rasters to remove outer NA valuse
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## combine forecasts together
forecast <- stack(forecast2050, forecast2099)

## plot raster forecasts
nam <- c(paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" "), paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" "))

theme <- rasterTheme(region=colorRampPalette(c("darkorchid4","darkorchid","darkorchid1","yellow","yellow2","gold","gold3"))(100)) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(4,6), # horizontal layout
               names.attr=nam, # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Food Safety Inspections Due to Future Warming\n(per 100,000 Facilities)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot
jpeg(filename = "./inspection_forecastRCP85.jpeg", width=8, height=8, units = "in", type="cairo", res=300)
p
dev.off()

# viol----
id2050 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set month id

id2099 <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set month id

forecast2050 <- setZ(diff2050.mn.viol*1000, id2050) # set year as z variable and convert to per 1k
names(forecast2050) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" ") # set names

forecast2099 <- setZ(diff2099.mn.viol*1000, id2099) # set year as z variable and convert to per 1k
names(forecast2099) <- paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" ") # set names

## trim forecast rasters to remove outer NA valuse
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

## combine forecasts together
forecast <- stack(forecast2050, forecast2099)

## plot raster forecasts
nam <- c(paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2050", sep=" "), paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December"), "2099", sep=" "))

theme <- rasterTheme(region=colorRampPalette(c("gray85","gray95","darkolivegreen1","darkolivegreen2","darkolivegreen3","darkolivegreen4","darkolivegreen"))(100)) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(4,6), # horizontal layout
               names.attr=nam, # inset titles
               colorkey=list(space="bottom", labels=list(cex=1)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Projected Change in Mean Daily Food Safety Violations Due to Future Warming\n(per 1,000 Inspections)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot
jpeg(filename = "./violation_forecastRCP85.jpeg", width=8, height=8, units = "in", type="cairo", res=300)
p
dev.off()
```

```{r pol.grids8.5, out.width="100%", fig.cap=paste0(sfig_num("pol.grids8.5", display="cite"),". Grid cell projections of change in number of police stops for RCP8.5. This figure presents the 25km x 25km grid cell projections of the potential impact of warming maximum temperatures on the percentage point change in police stops in future months. In this figure, downscaled climatic model maximum temperature projections are averaged across the 21 models in the ensemble and then coupled with our historical model parameters to produce an estimated percentage point change in police stops in each geographic location for the months of 2050 and 2099. Areas of the northern United States -- where non-summer temperatures are currently coldest -- may experience increases in net police regulatory activity due to future warming. However, southern portions of the U.S. may experience net reductions in police regulatory stops due to future warming.")}
knitr::include_graphics("./police_forecastRCP85.jpeg")
```

```{r cr.grids8.5, out.width="100%", fig.cap=paste0(sfig_num("cr.grids8.5", display="cite"),". Grid cell projections of change in number of fatal vehicular crashes for RCP8.5. This figure presents the 25km x 25km grid cell projections of the potential impact of warming maximum temperatures on the percentage point change in fatal crashes in future months. In this figure, downscaled climatic model maximum temperature projections are averaged across the 21 models in the ensemble and then coupled with our historical model parameters to produce an estimated percentage point change in fatal crashes in each geographic location for the months of 2050 and 2099. Areas of the northern United States -- where non-summer temperatures are currently coldest -- may experience decreases in net fatal crashes in the winter due to future warming. However, southern portions of the U.S. -- especially in summer months -- may experience net increases in fatal vehicular crashes due to future warming.")}
knitr::include_graphics("./crash_forecastRCP85.jpeg")
```

```{r insp.grids8.5, out.width="100%", fig.cap=paste0(sfig_num("insp.grids8.5", display="cite"),". Grid cell projections of change in food safety inspections for RCP8.5. This figure presents the 25km x 25km grid cell projections of the potential impact of warming maximum temperatures on the change in food safety regulatory inspections in future months. In this figure, downscaled climatic model maximum temperature projections are averaged across the 21 models in the ensemble and then coupled with our historical model parameters to produce an estimated change in food safety inspections per 1,000 facilities in each geographic location for the months of 2050 and 2099. Areas of the northern United States may experience increases in net food safety inspection activity during the winter months due to future warming. Southern portions of the U.S. may experience net reductions in inspection activity -- particularly during summer months when food safety risks are highest -- due to future warming.")}
knitr::include_graphics("./inspection_forecastRCP85.jpeg")
```

```{r viol.grids8.5, out.width="100%", fig.cap=paste0(sfig_num("viol.grids8.5", display="cite"),". Grid cell projections of change in number of food safety violations for RCP8.5. This figure presents the 25km x 25km grid cell projections of the potential impact of warming maximum temperatures on the change in food safety violations in future months. In this figure, downscaled climatic model maximum temperature projections are averaged across the 21 models in the ensemble and then coupled with our historical model parameters to produce an estimated change in food safety violations per 1,000 inspections for each geographic location for the months of 2050 and 2099. Our projection suggests that warming may increase food safety violation rates across the entire US by 2099.")}
knitr::include_graphics("./violation_forecastRCP85.jpeg")
```
